## ---------------------------
##
## Script name: Chicago_Analysis.R
##
## Description: Outputs all figures for Chicago Analysis.
##
## ---------------------------
## DEPENDENCIES

require(dplyr)
library(tidyverse)
library(ggplot2)
library(sf)
library(tigris)
library(stringr)
library(grid)
library(gtable)
library(kableExtra)
library(gridExtra)

#------------------- LOAD HELPER FUNCTIONS ------------------

source(paste(file_locations$current_fp, "/R/Chicago/Chicago_graphing_helpers.R", sep="/"))

#---------------------- LOAD DATA ---------------------------
# murder data
city_murder_data <- read.csv(paste(file_locations$current_fp, file_locations$Chicago$Murders$processed$city, sep="/"))

tract_murder_data <- read.csv(paste(file_locations$current_fp, file_locations$Chicago$Murders$processed$tract, sep="/")) %>%
    mutate(tract_2010 = str_pad(tract_2010, 11, side="left", pad="0"))

# demographic data
demos <- read_csv(paste(file_locations$current_fp, file_locations$Chicago$Demographics$processed, sep="/")) %>%
  mutate(tract_2010 = str_pad(tract_2010, 11, side="left", pad="0"))

#------------------------------------------------------------

# filepath to output
figure_fp <- paste(file_locations$current_fp, "figures/Chicago", sep="/")
options(knitr.table.format = "latex")

#------------------------------------------------------------

run_Chicago_analysis <- function(file_locations) {

  #---------------------------------------------
  # 1. OVERALL CHICAGO MURDER RATES
  #---------------------------------------------
  
  ggplot(city_murder_data, aes(x=year, y=murder_rate)) +
    geom_line() +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    ylim(0, 40) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 100K residents") +
    # ggtitle("MURDER RATE IN CHICAGO\
    # 1965-2020") +
    theme(axis.title.x = element_text(size = 16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20))
  
  ggsave(paste(figure_fp, "Figure1_overall_murders_1965-2020_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  #---------------------------------------------
  # 2. Chicago Murder Rates - 1965-2020 - Maps
  #---------------------------------------------

  ### A. Plot Maps
  tracts_df <- tracts(state = "IL", year=2010, cb=TRUE)
  tracts_df <- st_as_sf(tracts_df) %>% st_transform(crs=4326)
  
  tract_murders_long <- tract_murder_data %>%
    dplyr::select(city, tract_2010, starts_with("murders") & !contains("2014.2020")) %>%
    gather(key=year, value=total_murders, starts_with("murders")) %>%
    mutate(year = str_split(year, "_", simplify=T)[,2],
           year = str_replace_all(year, "[.]", "-"))
  
  tract_murders_long <- tracts_df %>%
    unite("tractname", STATE, COUNTY, TRACT, sep="") %>%
    inner_join(tract_murders_long, by=c("tractname"="tract_2010"))
  
  # split data for each "facet" i.e. time period
  df <- tract_murders_long %>%
    filter(year %in% c("1965-1970","1971-1975","1976-1980","1981-1985"))

  df <- split(df,f = df$year)
  
  p1 <- ggplot(df[[1]]) +
    geom_sf(aes(fill=total_murders), size=0.01,  color="gray75") +
    facet_wrap(~year, ncol = 1, labeller = labeller(VAR = "year")) +
    theme_minimal() +
    theme(axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major = element_blank(),
          legend.position="bottom",
          legend.title = element_blank(),
          legend.key.size = unit(0.4, 'cm'),
          legend.text = element_text(size=8),
          strip.text = element_text(size=10)) +
    scale_fill_gradient(low = "white", high = "indianred4", limits=c(0,62))
  
  # do it for each "facet"
  p2 <- p1 %+% df[[2]]
  p3 = p2 %+% df[[3]]
  p4 = p3 %+% df[[4]]
  
  # combine facets
  grobs_homs1 <- (arrangeGrob(p1 + theme(legend.position="none"),
                                p2 + theme(legend.position="none"),
                                p3 + theme(legend.position="none"),
                                p4 + theme(legend.position="none"), nrow=1))
  
  # do same thing for time periods in lower portion of the figure
  df <- tract_murders_long %>%
    filter(year %in% c("2001-2005","2006-2010","2011-2015","2016-2020"))
  
  df <- split(df,f = df$year)
  
  p1 <- ggplot(df[[1]]) +
    geom_sf(aes(fill=total_murders), size=0.01, color="gray75") +
    facet_wrap(~year, ncol = 1, labeller = labeller(VAR = "year")) +
    theme_minimal() +
    theme(axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major = element_blank(),
          legend.title = element_blank(),
          legend.key.size = unit(0.4, 'cm'),
          legend.text = element_text(size=8),
          strip.text = element_text(size=10)) +
    scale_fill_gradient(low = "white", high = "indianred4", limits=c(0,62))
  
  # do it for each "facet"
  p2 <- p1 %+% df[[2]]
  p3 = p2 %+% df[[3]]
  p4 = p3 %+% df[[4]]
  
  # combine facets
  grobs_homs2 <- (arrangeGrob(p1 + theme(legend.position="none"),
                                p2 + theme(legend.position="none"),
                                p3 + theme(legend.position="none"),
                                p4 + theme(legend.position="none"), nrow=1))
  
  # create a single uniform legend for all maps
  g_legend <- function(a.gplot){
    tmp <- ggplot_gtable(ggplot_build(a.gplot))
    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
    legend <- tmp$grobs[[leg]]
    return(legend)}
  
  mylegend <- g_legend(p1)
  
  # set output file for the plot
  pdf(paste(figure_fp, "Figure2_chicago_murder_maps.pdf", sep="/"), 
      width = 7, height = 4)
  
  # arrange and draw the figure for output
  grid.draw(grid.arrange(arrangeGrob(grobs_homs1,
                           grobs_homs2,
                           nrow=2),
               mylegend, nrow=1, ncol=2, heights=c(1), widths=c(10, 2)))
  
  dev.off()
  
  #---------------------------------------------
  # 3. Chicago Murder Rates by 1965-1970 Quintile
  #---------------------------------------------
  city_tracts <-  read_sf(paste(file_locations$current_fp, file_locations$Nation$Geographies$processed, sep="/")) %>%
    filter(city == "Chicago") 
  
  # many tracts have zero murders but will still be ranked
  # and assigned a quintile; so we will just arrange the 
  # tracts numerically
  tract_murder_data_ordered <- city_tracts %>%
    st_drop_geometry(.) %>%
    dplyr::select(city, tract_2010) %>%
    inner_join(tract_murder_data, by=c("city","tract_2010"))
  
  # determine 1965-1970 murder rate quintile
  quintile_labels <- tract_murder_data_ordered %>%
    dplyr::select(city, tract_2010, murder_rate_1965.1970) %>%
    arrange(murder_rate_1965.1970) %>%
    mutate(nrows = nrow(.),
           t = as.numeric(rownames(.)),
           quintile = case_when(t >= 1 & t <= floor(nrows/5) ~ "First",
                                t > floor(nrows/5) & t <= floor((nrows/5)*2) ~ "Second",
                                t > floor((nrows/5)*2) & t <= floor((nrows/5)*3) ~ "Third",
                                t > floor((nrows/5)*3) & t <= floor((nrows/5)*4) ~ "Fourth",
                                t > floor((nrows/5)*4)  ~ "Fifth")) %>%
    select(tract_2010, quintile) %>% distinct()
  
  # gather all 5-year increment murder rates
  tract_murders_long <- tract_murder_data_ordered %>%
    dplyr::select(city, tract_2010, starts_with("murder_rate") & !contains("2014.2020")) %>%
    gather(key=year, value=murder_rate, starts_with("murder_rate")) %>%
    mutate(year = str_split(year, "rate_", simplify=T)[,2],
           year = str_split(year, "[.]", simplify=T)[,2],
           year = as.numeric(year))

    
  subdf <- tract_murders_long %>%
    inner_join(quintile_labels, by="tract_2010") %>%
    group_by(quintile, year) %>%
    summarize(avg_murder_rate = mean(murder_rate, na.rm=T))
  
  subdf$quintile <- factor(subdf$quintile, levels=c("First","Second","Third","Fourth","Fifth"))

  ggplot(subdf, aes(x=year, y=avg_murder_rate, linetype=quintile)) +
    geom_line(aes(group=quintile)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Average Murder rate per 10K residents") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.8, 0.8),
          legend.text = element_text(size = 16))
  
  ggsave(paste(figure_fp, "Figure3_homicide_avg_nbhd_quintiles_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  
  #----------------------------------------------------------
  # 4. Neighborhood Murder Rates by Poverty Level & Race
  #----------------------------------------------------------
  
  # 4a. BY POVERTY
  avg_rates <- city_murder_data %>%
    dplyr::select(city, year, avg_poor_rate, avg_nonpoor_rate) %>%
    gather(key=poverty_bin, value=avg_murder_rate, avg_poor_rate, avg_nonpoor_rate) %>%
    mutate(poverty_bin = str_split(poverty_bin, "_", simplify=T)[,2])
  
  avg_rates$poverty_bin <- factor(avg_rates$poverty_bin, levels=c("poor","nonpoor"))
  
  ggplot(avg_rates, 
         aes(x=year, y=avg_murder_rate, linetype=poverty_bin)) +
    geom_line() +
    ylim(c(0,10)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 10K residents") +
    # ggtitle("NEIGHBORHOOD VIOLENCE IN MAJORITY-BLACK AND\
    # MAJORITY-WHITE NEIGHBORHOODS\
    # 1965-2020") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.55, 0.93),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Average murder rate in poor neighborhoods",
                                   "Average murder rate in nonpoor neighborhoods"))
  
  ggsave(paste(figure_fp, "Figure4a_homicide_avg_nbhd_poverty_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  
  # Diff/Ratios
  avg_rates_2 <- avg_rates %>%
    spread(key = poverty_bin, value=avg_murder_rate) %>%
    mutate(pnp_diff = poor - nonpoor,
           pnp_ratio = poor/nonpoor) %>%
    gather(key="var", value="value", pnp_diff, pnp_ratio)
  
  avg_rates_2$var <- factor(avg_rates_2$var, levels=c("pnp_diff","pnp_ratio"))
  
  ggplot(avg_rates_2, aes(x=year, y=value, linetype=var)) +
    geom_line() +
    ylim(c(0,10)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murder rate per 10K residents \n in poor vs. nonpoor neighborhoods") +
    # ggtitle("NEIGHBORHOOD VIOLENCE IN MAJORITY-BLACK AND\
    # MAJORITY-WHITE NEIGHBORHOODS\
    # 1965-2020") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.5, 0.87),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Difference in murder rate in poor vs. \n nonpoor neighborhoods",
                                   "Ratio of murder rate in poor vs. \n nonpoor neighborhoods"))
  
  ggsave(paste(figure_fp, "Figure4b_homicide_avg_nbhd_poverty_Chicago_ratiodiff.png", sep="/"),
         width = 8, height = 6)
  
  # 4b. BY RACE
  
  avg_rates <- city_murder_data %>%
    dplyr::select(city, year, avg_black_rate, avg_white_rate) %>%
    gather(key=majority_race_bin, value=avg_murder_rate, avg_black_rate, avg_white_rate) %>%
    mutate(majority_race_bin = str_split(majority_race_bin, "_", simplify=T)[,2])

  ggplot(avg_rates, aes(x=year, y=avg_murder_rate, linetype=majority_race_bin)) +
    geom_line() +
    ylim(c(0,10)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 10K residents") +
    # ggtitle("NEIGHBORHOOD VIOLENCE IN MAJORITY-BLACK AND\
    # MAJORITY-WHITE NEIGHBORHOODS\
    # 1965-2020") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.45, 0.85),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Average murder rate in majority-Black neighborhoods",
                                   "Average murder rate in majority-White neighborhoods"))
  
  ggsave(paste(figure_fp, "Figure5a_homicide_avg_nbhd_race_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  # Diff/Ratios
  
  avg_rates_2 <- avg_rates %>%
    select(year, majority_race_bin, avg_murder_rate) %>%
    spread(key = majority_race_bin, value=avg_murder_rate) %>%
    mutate(bw_diff = black - white,
           bw_ratio = black/white) %>%
    gather(key="var", value="value", bw_diff, bw_ratio)
  
  avg_rates_2$var <- factor(avg_rates_2$var, levels=c("bw_diff","bw_ratio"))
  
  ggplot(avg_rates_2, aes(x=year, y=value, linetype=var)) +
    geom_line() +
    scale_y_continuous(breaks = seq(0, 10, by = 1)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murder rate per 10K residents in\nmajority-Black vs. majority-White neighborhoods") +
    # ggtitle("NEIGHBORHOOD VIOLENCE IN MAJORITY-BLACK AND\
    # MAJORITY-WHITE NEIGHBORHOODS\
    # 1965-2020") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.37, 0.87),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Difference in murder rate in majority-Black vs. \n majority-White neighborhoods",
                                   "Ratio of murder rate in majority-Black vs. \n majority-White neighborhoods"))
  
  ggsave(paste(figure_fp, "Figure5b_homicide_avg_nbhd_race_Chicago_ratiodiff.png", sep="/"),
         width = 8, height = 6)
  
  #----------------------------------------------------------
  # 5. Exposure Rates by Poverty Level & Race
  #----------------------------------------------------------
  
  # 5a. BY POVERTY
  
  avg_exposures <- city_murder_data %>%
    dplyr::select(city, year, avg_poor_exp_rate, avg_nonpoor_exp_rate) %>%
    gather(key=poverty_bin, value=avg_murder_exp_rate, avg_poor_exp_rate, avg_nonpoor_exp_rate) %>%
    mutate(poverty_bin = str_split(poverty_bin, "_", simplify=T)[,2])
  
  avg_exposures$poverty_bin <- factor(avg_exposures$poverty_bin, levels=c("poor", "nonpoor"))
  
  ggplot(avg_exposures, aes(x=year, y=avg_murder_exp_rate, linetype=poverty_bin)) +
    geom_line() +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    ylim(c(0,10)) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 10K residents") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.5, 0.8),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Murder rate in neighborhood of average poor resident",
                                   "Murder rate in neighborhood of average nonpoor resident"))
  
  ggsave(paste(figure_fp, "Figure6a_homicide_exposure_poverty_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  # Diff/Ratios
  
  avg_exposures_2 <- avg_exposures %>%
    spread(key = poverty_bin, value=avg_murder_exp_rate) %>%
    mutate(diff = poor - nonpoor,
           ratio = poor/nonpoor) 
  
  avg_exposures_2 <- avg_exposures_2 %>%
    gather(key="var", value="value", diff, ratio)
  
  ggplot(avg_exposures_2, aes(x=year, y=value, linetype=var)) +
    geom_line() +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    ylim(c(0,5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murder rate per 10K residents \n in neighborhoods of poor vs nonpoor residents") +
    theme(axis.title.x = element_text(margin = margin(t = 10, r = 20, b = 0, l = 0),
                                      size=16, face="bold"),
          axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0),
                                      size=16, face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.5, 0.8),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Difference in murder rate in neighborhood of average \n poor vs nonpoor resident",
                                   "Ratio of murder rate in neighborhood of average poor \n vs nonpoor resident"))
  
  ggsave(paste(figure_fp, "Figure6b_homicide_exposure_poverty_Chicago_ratiodiff.png", sep="/"),
         width = 8, height = 6)
  
  # 5b. BY RACE
  
  avg_exposures <- city_murder_data %>%
    dplyr::select(city, year, avg_black_exp_rate, avg_white_exp_rate, avg_otherrace_exp_rate) %>%
    gather(key=race_bin, value=avg_murder_exp_rate, avg_black_exp_rate, avg_white_exp_rate, avg_otherrace_exp_rate) %>%
    mutate(race_bin = str_split(race_bin, "_", simplify=T)[,2])

  ggplot(avg_exposures, aes(x=year, y=avg_murder_exp_rate, linetype=race_bin)) +
    geom_line() +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    ylim(c(0,10)) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 10K residents") +
    # ggtitle("EXPOSURE TO NEIGHBORHOOD VIOLENCE\
    # FOR BLACK AND WHITE RESIDENTS\
    # 1965-2020") +
    theme(axis.title.x = element_text(margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold", size=16),
          axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold", size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          legend.title = element_blank(),
          legend.position = c(0.5, 0.87),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed", "longdash"),
                          labels=c("Murder rate in neighborhood of average black resident",
                                   "Murder rate in neighborhoods of average residents \n from all other ethnic groups",
                                   "Murder rate in neighborhood of average white resident"))
  
  ggsave(paste(figure_fp, "Figure7a_homicide_exposure_race_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  # Diff/Ratios
  
  avg_exposures_2 <- avg_exposures %>%
    spread(key = race_bin, value=avg_murder_exp_rate) %>%
    mutate(bw_diff = black - white,
           bw_ratio = black/white) 
  
  avg_exposures_2 <- avg_exposures_2 %>%
    gather(key="var", value="value", bw_diff, bw_ratio)
  
  ggplot(avg_exposures_2, aes(x=year, y=value, linetype=var)) +
    geom_line() +
    ylim(c(0,10)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murder rate per 10K residents in\nneighborhoods of average Black vs. White\nresident") +
    # ggtitle("EXPOSURE TO NEIGHBORHOOD VIOLENCE\
    # FOR BLACK AND WHITE RESIDENTS\
    # 1965-2020") +
    theme(axis.title.x = element_text(margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold", size=16),
          axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold", size=16),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.5, 0.87),
          legend.text = element_text(size=16)) +
    scale_linetype_manual(values=c("solid", "dashed"),
                          labels=c("Difference in murder rate in neighborhood of average \n Black vs White resident",
                                   "Ratio of murder rate in neighborhood of average \n Black vs White resident"))
  
  ggsave(paste(figure_fp, "Figure7b_homicide_exposure_race_Chicago_ratiodiff.png", sep="/"),
         width = 8, height = 6)
  
  #----------------------------------------------------------
  # 6. Exposure Rates for Black Residents In/Outside Majority-Black neighborhoods
  #----------------------------------------------------------
  
  avg_exposures <- city_murder_data %>%
    dplyr::select(city, year, avg_black_exp_rate, avg_majority_black_exp_rate, 
                  avg_nonmajority_black_exp_rate) %>%
    gather(key=majority_race_bin, value=avg_murder_exp_rate, avg_black_exp_rate, avg_majority_black_exp_rate, 
           avg_nonmajority_black_exp_rate) %>%
    mutate(majority_race_bin = str_split(majority_race_bin, "_", simplify=T)[,2])
  
  ggplot(avg_exposures, aes(x=year, y=avg_murder_exp_rate, linetype=majority_race_bin)) +
    geom_line() +
    ylim(c(0,10)) +
    scale_x_continuous(breaks = seq(1965, 2020, by = 5)) +
    theme_light() +
    xlab("Year") +
    ylab("Murders per 10K residents") +
    # ggtitle("EXPOSURE TO NEIGHBORHOOD VIOLENCE FOR\
    # BLACK RESIDENTS IN AND OUTSIDE OF BLACK\
    # NEIGHBORHOODS, 1965-2020") +
    theme(axis.title.x = element_text(size=16, margin = margin(t = 10, r = 20, b = 0, l = 0), face="bold"),
          axis.title.y = element_text(size=16, margin = margin(t = 0, r = 20, b = 0, l = 0), face="bold"),
          axis.text.x = element_text(size=16),
          axis.text.y = element_text(size=16),
          plot.title = element_text(face="bold", hjust = 0.5, size=20),
          legend.title = element_blank(),
          legend.position = c(0.45, 0.85),
          legend.text = element_text(size = 16)) +
    scale_linetype_manual(values=c("solid", "dashed", "longdash"),
                          labels=c("Murder rate in neighborhoods of average Black resident",
                                   "Murder rate in neighborhoods of average Black resident \n in majority-Black neighborhood",
                                   "Murder rate in neighborhoods of average Black resident \n NOT in majority-Black neighborhood"))
  
  
  ggsave(paste(figure_fp, "Figure8_homicide_exposure_black_Chicago.png", sep="/"),
         width = 8, height = 6)
  
  
  #----------------------------------------------------------
  # 7. Compounded Disadvantage
  #----------------------------------------------------------

  # 7a. Murders, Incarceration, Police Shootings MAPS
  
  # MURDERS
  murder_rates <- tract_murder_data %>%
    dplyr::select(city, tract_2010, murder_rate_2014.2020)
  
  tracts_df <- tracts(state = "IL", year=2010, cb=TRUE)
  tracts_df <- st_as_sf(tracts_df) %>% st_transform(crs=4326)
  
  agg_df <- tracts_df %>%
    unite("tract_2010", STATE, COUNTY, TRACT, sep="") %>%
    inner_join(murder_rates, by=c("tract_2010")) %>%
    mutate(year = "2014-2020")
  
  # split data for each "facet"
  df_plot <- split(agg_df,f = agg_df$year)
  
  # plot for the first "facet"
  p_homs <- ggplot(df_plot[[1]]) +
    geom_sf(aes(fill=murder_rate_2014.2020), size=0.1, color="gray75") +
    facet_wrap(~year, ncol = 1, labeller = labeller(VAR = "year")) +
    theme_minimal() +
    ggtitle("Murder Rate") +
    theme(axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major = element_blank(),
          legend.position="bottom",
          legend.title = element_blank(),
          legend.key.size = unit(0.6, 'cm'),
          legend.text = element_text(size=8),
          strip.text = element_text(size=10),
          plot.title = element_text(size=11, hjust=0.5))+
    scale_fill_gradient(low = "white", high = "indianred4")
  
  
  # POLICE SHOOTINGS
  ps_df <- read.csv(paste(file_locations$current_fp, file_locations$Chicago$`Police Shootings`$processed, sep="/")) %>%
    mutate(tract_2010 = str_pad(tract_2010, 11, side="left", pad="0")) %>%
    filter(city=="Chicago")
  
  ps_df <- demos %>% filter(year %in% ps_df$year) %>%
    inner_join(ps_df, by=c("tract_2010", "year"))

  ps_rates <- .calc_crime_rate(df=ps_df, per=10000, crime_col="fatal_police_shootings", population_col="total_population", 
                   grouping_vars=c("city","tract_2010")) %>%
    rename(ps_rate = crime_rate) %>%
    mutate(year = "2014-2020")
  
  tracts_df <- tracts(state = "IL", year=2010, cb=TRUE)
  tracts_df <- st_as_sf(tracts_df) %>% st_transform(crs=4326)
  
  agg_df <- tracts_df %>%
    unite("tract_2010", STATE, COUNTY, TRACT, sep="") %>%
    inner_join(ps_rates, by=c("tract_2010"))
  
  # split data for each "facet"
  df_plot <- split(agg_df,f = agg_df$year)
  
  # plot for the first "facet"
  p_ps <- ggplot(df_plot[[1]]) +
    geom_sf(aes(fill=ps_rate), size=0.1, color="gray75") +
    facet_wrap(~year, ncol = 1, labeller = labeller(VAR = "year")) +
    theme_minimal() +
    ggtitle("Police Shooting Rate") +
    theme(axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major = element_blank(),
          legend.position="bottom",
          legend.title = element_blank(),
          legend.key.size = unit(0.6, 'cm'),
          legend.text = element_text(size=8),
          strip.text = element_text(size=10),
          plot.title = element_text(size=11, hjust=0.5))+
    scale_fill_gradient(low = "white", high = "blue4")
  
  # MALE INCARCERATION
  inc_rates <- read.csv(paste(file_locations$current_fp, file_locations$Chicago$`Male Incarceration`$processed, sep="/")) %>%
    mutate(tract_2010 = str_pad(tract_2010, 11, side="left", pad="0"))
  
  inc_rates <- demos %>% filter(year %in% inc_rates$year) %>%
    left_join(inc_rates, by=c("tract_2010", "year")) %>%
    mutate(incarc_rate = ifelse(is.na(incarc_rate), 0, incarc_rate),
           year=2010)
  
  tracts_df <- tracts(state = "IL", year=2010, cb=TRUE)
  tracts_df <- st_as_sf(tracts_df) %>% st_transform(crs=4326)
  
  agg_df <- tracts_df %>%
    unite("tract_2010", STATE, COUNTY, TRACT, sep="") %>%
    inner_join(inc_rates, by=c("tract_2010"))
  
  # split data for each "facet"
  df_plot <- split(agg_df,f = agg_df$year)
  
  p_incarc <- ggplot(df_plot[[1]]) +
    geom_sf(aes(fill=incarc_rate), size=0.1, color="gray75") +
    facet_wrap(~year, ncol = 1, labeller = labeller(VAR = "year")) +
    theme_minimal() +
    ggtitle("Male Incarceration Rate") +
    theme(axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          panel.grid.major = element_blank(),
          legend.position="bottom",
          legend.title = element_blank(),
          legend.key.size = unit(0.6, 'cm'),
          legend.text = element_text(size=8),
          strip.text = element_text(size=10),
          plot.title = element_text(size=11, hjust=0.5)) +
    scale_fill_gradient(low = "white", high = "springgreen4",
                        breaks = seq(0, 2000, 1000))
  
  g <- arrangeGrob(p_homs, p_incarc, p_ps, nrow=1)
  
  pdf(paste(figure_fp, "Figure9_homs_ps_incarc_spatial_plots.pdf", sep="/"),
      width=7, height=5)
  
  grid.draw(g)
  
  dev.off()
  
  
  # 7b. Compounded Disadvantage by Groups
  labels1 <- .generate_tract_labels_Chicago(demos, 2014, 2020)
  murder_quintile_labels1 <- tract_murder_data %>%
    dplyr::select(city, tract_2010, murder_rate_2014.2020) %>%
    arrange(murder_rate_2014.2020) %>%
    mutate(nrows = nrow(.),
           breaks <- floor(nrows / 5),
           t = as.numeric(rownames(.)),
           quintile = case_when(t >= 1 & t <= breaks ~ "First",
                                t > breaks & t <= 2*breaks ~ "Second",
                                t > 2*breaks & t <= 3*breaks ~ "Third",
                                t > 3*breaks & t <= 4*breaks ~ "Fourth",
                                t > 4*breaks ~ "Fifth")) %>%
    mutate(quintile = ifelse(quintile%in%c("Second","Third","Fourth"), "middle_quintiles", 
                             paste0(str_to_lower(quintile),"_quintile"))) %>%
    select(tract_2010, quintile) %>% distinct()

  labels2 <- .generate_tract_labels_Chicago(demos, 1965, 1970)
  murder_quintile_labels2 <- tract_murder_data %>%
    dplyr::select(city, tract_2010, murder_rate_1965.1970) %>%
    arrange(murder_rate_1965.1970) %>%
    mutate(nrows = nrow(.),
           breaks <- floor(nrows / 5),
           t = as.numeric(rownames(.)),
           quintile = case_when(t >= 1 & t <= breaks ~ "First",
                                t > breaks & t <= 2*breaks ~ "Second",
                                t > 2*breaks & t <= 3*breaks ~ "Third",
                                t > 3*breaks & t <= 4*breaks ~ "Fourth",
                                t > 4*breaks ~ "Fifth")) %>%
    mutate(quintile = ifelse(quintile%in%c("Second","Third","Fourth"), "middle_quintiles", 
                             paste0(str_to_lower(quintile),"_quintile"))) %>%
    select(tract_2010, t, quintile) %>% distinct()
  
  tract_labels <- (labels1 %>% inner_join(murder_quintile_labels1, by="tract_2010")) %>%
    inner_join((labels2 %>% inner_join(murder_quintile_labels2, by="tract_2010")), 
               by=c("tract_2010"), suffix=c("_14_20", "_65_70")) %>%
    mutate(city="Chicago")
  
  # 2010-2020 TABLE
  
  df_inc_14_20 <- .generate_tables_by_groups_Chicago(rates_df=inc_rates, rate_col="incarc_rate", label_suffix="_14_20",
                                             tract_labels) %>%
    mutate(incarc_rate = round(incarc_rate, 0))
  df_ps_14_20 <- .generate_tables_by_groups_Chicago(rates_df=ps_rates, rate_col="ps_rate", label_suffix="_14_20",
                                            tract_labels)
  
  subdf_14_20 <- df_inc_14_20 %>%
    inner_join(df_ps_14_20, by=c("category","group")) %>%
    arrange(category, group)
  
  subdf_14_20$category <- plyr::revalue(subdf_14_20$category, c("Quintile"="2014-2020 Quintile",
                                                                "Race"="2014-2020 Race",
                                                                "Poverty"="2014-2020 Poverty"))
  
  landscape(kable(subdf_14_20%>%select(-category), 
      align=c("l","c", "c"), longtable=T, booktabs=TRUE,
      col.names = c( "", "2010 Male Incarceration Rate","2014-2020 Police Shootings Rate")) %>%
    column_spec(1, italic = T, width="3cm") %>%
    row_spec(0, bold=T) %>%
    column_spec(2:3, width="7cm") %>%
    pack_rows(index = table(subdf_14_20$category)) %>%
    kable_classic(full_width = F,
                  font_size = 10, latex_options = c("striped", "repeat_header"),
                  position = "center")) %>%
    save_kable(paste(figure_fp, "Table1_comp_disadv_table_2014-2020Quintiles.pdf", sep="/"))
  
  
  # 1965-1970 TABLE
  
  df_inc_65_70 <- .generate_tables_by_groups_Chicago(rates_df=inc_rates, rate_col="incarc_rate", label_suffix="_65_70",
                                             tract_labels) %>%
    mutate(incarc_rate = round(incarc_rate, 0))
  df_ps_65_70 <- .generate_tables_by_groups_Chicago(rates_df=ps_rates, rate_col="ps_rate", label_suffix="_65_70",
                                            tract_labels)
  
  subdf_65_70 <- df_inc_65_70 %>%
    inner_join(df_ps_65_70, by=c("category","group")) %>%
    arrange(category, group)
  
  subdf_65_70$category <- plyr::revalue(subdf_65_70$category, c("Quintile"="1965-1970 Quintile",
                                                                "Race"="1965-1970 Race",
                                                                "Poverty"="1965-1970 Poverty"))

  landscape(kable(subdf_65_70%>%select(-category), 
      align=c("l","c", "c"), longtable=T, booktabs=TRUE,
      col.names = c( "", "2010 Male Incarceration Rate","2014-2020 Police Shootings Rate")) %>%
    column_spec(1, italic = T, width="3cm") %>%
    row_spec(0, bold=T) %>%
    column_spec(2:3, width="7cm") %>%
    pack_rows(index = table(subdf_65_70$category)) %>%
    kable_classic(full_width = F,
                  font_size = 10, latex_options = c("striped", "repeat_header"),
                  position = "center")) %>%
    save_kable(paste(figure_fp, "Table2_comp_disadv_table_1965-1970Quintiles.pdf", sep="/"))
}

